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PREFACE 


I  prepared  this  project  as  an  article  for  publication.  I 
wrote  it  in  the  style  of  the  target  journal  which  explains  the 
use  of  double  spacing,  the  heavy  use  of  passive  voice,  and  the 
unusual  format  for  the  references.  Subject  to  clearance,  this 
manuscript  will  be  submitted  to  The  International  Journal  for 
Numerical  Methods  in  Engineering  for  consideration. 

There  are  parallel  methods  for  solving  parabolic 
differential  equations  in  both  finite  elements  and  finite 
differences.  Quite  often,  the  method  of  choice  for  treating  the 
temporal  domain  is  the  Crank-Nicolson  method.  A  closely  related 
method  for  finite  differences  is  the  Crandall  method. 

Similarities  between  the  matrix  equations  for  the  Crank-Nicolson 
versions  of  finite  differences  and  finite  elements  suggested 
searching  for  an  equivalent  for  the  Crandall  method  for  finite 
elements.  I  found  such  a  method  to  exist,  and  it  is  derived 
herein.  I  have  also  provided  a  theoretical  treatment  of  accuracy 
and  stability  along  with  numerical  results  for  validation. 
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AN  OPTIMUM  FORMULATION  OF  THE  FINITE  ELEMENT 
METHOD  FOR  THE  DIFFUSION  EQUATION 


C.  R.  MARTIN 

Air  Command  and  Staff  College,  Maxwell  AFB,  Alabama,  U.S.A. 


SUMMARY 

An  optimum  formulation  of  the  finite  element  method  for  the 
diffusion  equation  is  derived  and  applied  to  a  problem  with 
Dirichlet  and  Neumann  boundary  conditions.  The  method  is  the 
finite  element  analog  of  the  Crandall  method1  for  finite 
differences . 


INTRODUCTION 


In  an  earlier  paper,  the  Crandall  HOT!  and  Crartk-Nicolson 
(CNM)  finite  difference  methods  were  compared  with  respect  to 
the  solution  of  the  transient  heat  conduction  problem  (isomorphic 
with  respect  to  the  diffusion  equation). 

In  this  paper,  an  optimum  formulation  of  the  finite  element 
method  FEM  h  for  the  diffusion  equation  is  derived.  This  method 
is  the  finite  element  equivalent  of  the  Crandall  method  for 
finite  differences.  The  method  is  then  applied  to  problems  with 
Neumann  and  Dirichlet  boundary  conditions,  and  the  accuracy  is 
compared  with  the  finite  element  version  of  the  CNM. 
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DERIVATION  OF  THE  METHOD 


The  following  quadratic  functional,  given  by  Myers, 

I(,f)  =  ^0[(f£)2  +  T*  ♦*]  dx 

is  related  to  the  normalized  diffusion  equation 

2£  . 

38  1? 


with  boundary  conditions 


3$ 

'5x 


<K0,  e)  =  o  ,  e  >  6 

-  0  ,  0  >  0 


x=l 


and  initial  condition 


where 


0(xf  0)  *  1 


(1) 


(2) 


$  is  the  normalized  flux 
0  is  the  normalized  time  variable 
x  is  the  normalized  space  variable 

by  the  fact  that 

I(*)  i  I($)  for  all  f  in  H*  (3) 

where  0  is  the  solution  to  equation  (2)  and  where  Mg ,  the 
solution  space  for  equation  (2),  is  a  Hilbert  space  in  which  the 
member  functions  and  their  first  derivatives  are  square 
integrable  and  which  satisfy  the  essential  boundary  condition  of 
equation  (2). 

In  practice,  a  solution  is  sought  within  a  finite 
dimensional  subspace  of  Mg.  This  solution  will  be  the  Ritz 


E 

projection  of  the  true  soLution  onto  the  subspace  S  ,  where  the 
superscript  E  indicates  the  number  of  elements  used  in  the  finite 
element  solution. 

The  elements  of  s  are  the  trial  functions,  <j>  .  Since  they 
belong  to  Hg,  they  satisfy 

*E(0,  6)  =  0  ,  0  >  0  (4) 

The  simplest  choice  for  the  functions  <p  are  functions  which  are 

linear  within  each  element,  continuous  at  the  elemental 

_E 

interfaces,  and  zero  at  x  =  0.  Thus  $  is  expanded  in  terms  of 
the  familiar  tent  functions  (0  held  fixed): 

g 

iE(x)  =  *  T  S  rE(x)  (s) 

1  e=l  e  e 
E 

where  Te(*)  varies  linearly  from  (e-1)  Ax  to  eAx,  is  equal  to 
unity  at  eAx,  varies  linearly  from  eAx  to  (e  +  l)Ax,  and  is  equal 
to  zero  at  all  other  points  on  the  domain  including  (e-1)  Ax  and 
(e+l)Ax  •  The  interval  spacing  is  taken  to  be  constant,  thus 
Ax  =  1/E. 

The  problem,  then,  is  to  minimize 

'«*>  ‘  ^[(l?)2  *  ^V]dx  (6) 


over  a  single  element,  then 
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is  the  element  mass  matrix.  For  the  whole 
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where  0  =  •••  ?E1 >  K  is  the  global  stiffness  matrix,  and 

M  is  the  global  mass  matrix.  The  minimum  occurs  at  the  point 

T 

9l/3^j  *  0  for  the  vector  £  *  [(j>^  $2  •••  <J>g] 

determined  by 


*  =  0 


(9) 


If  equation  (9)  is  approximated  by  the  Crank-Nicolson  scheme, 
the  result  is 

/<Lk+1  -  (<Lk+1  *  <tk  \ 

T~  )-°  (10) 

An  optimum  formulation  can  be  derived  if  a  more  general  form 
of  equation  (10)  is  assumed 

(M  ♦  KaAt)  $_k+1  -  (M  -  K(l-a)  At)  <tk  C11) 


where  the  parameter  a  represents  the  weight  placed  on  the  next 
time  level,  k+1.  While  the  L*,  norm  is  not  a  natural  one  for  the 
finite  element  method,  it  is  of  interest  to  investigate  the  error 
with  respect  to  that  norm.  If  the  set  of  equations  represented 
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by  equation  (11)  are  treated  as  difference  equations,  then  the 

inherent  truncation  error  can  be  investigated  by  representing  the 

exact  solution  to  4  =  <ba  at  adjacent  points  in  the  space  time 

xx  o 

network  by  means  of  Taylor's  series  centered  at  (iAx,  kA0)  .  When 
these  are  substituted  into  the  general  expression 


A1  *i- 1 ,k+l  +  A2  *i,k+l  +  A1  *i+l,k+] 


(12) 


=  B1  *i-l.k  +  B2  ♦i.k  +  B1  *1*1.1 


where 


Aj  =  l-6ap 

A2  =  4+12ap 
B^  =  l+6(l-a)p 
B2  =  4-12 (1 -a) p 
p  =  A0/(Ax)2 

(p  is  known  within  the  literature  as  the  Fourier  modulus)  then 
the  remainder  terms  of  the  Taylor's  series  are 


Ax 


:(f  '  *  nr) 


yr ;  <J>  +  Ax*  p2(l-3a)  + 

Az;y  80 


(p(i  -  i)  +  rs) 


000 


♦  0  (Ax°)  + 
where  the  identities 


(13) 


*60  "  *exx  =  *xxxx  (14) 

090  *  0  0xx  *0xxxx  *xxxxxx 
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have  been  used.  It  is  obvious  that  the  expression  (11)  is  0(Ax  ) 
accurate  at  the  nodes.  An  0(Ax*)  scheme  can  be  derived  by 

7 

setting  the  coefficient  of  the  0(Ax  )  term  to  zero.  Thus  when 
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is  satisfied,  equation  (1L)  is  0(£x  j  accurate  at  the  nodes. 

2  4 

Furthermore,  at  the  intersection  point  of  the  0(&xL)  ana  0(Ax  j 
terms,  that  is  at  the  point 
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6 

the  truncation  error  of  (13)  and  therefore  of  (11)  is  O(AX^). 

The  stability  of  the  method  can  be  easily  investigated  for  a 
pure  Dirichlet  probLem  by  noting  that  equation  (11)  can  be 
written  as 


A  $_k  +  1  =  B  £k  (16) 

In  this  formulation,  the  first  and  last  equations  in  the  system 
nay  be  dropped,  thus  eliminating  the  known  boundary  values  and 
reducing  the  number  of  unknowns  by  two.  These  changes  will  not 
alter  the  numerical  solution.  In  the  resulting  system,  A.  and  B. 

are  tridiagonal.  Smi  th15  gives'  the  eigenvalues  for  the  iteration 

-  ^ 

matrix  C  =  A  B  as 


(l+2pa)  +  2(-pa)  cos  — 

(x  1  =  ..  _ N+l 

^c  n  ( 1  / . 

(l  +  2p(la))  *  2  (p  (1  -a ) )  cos  ^ 


where  N  is  the  number  of  equations.  The  analysis  will  be  done 
for  the  limiting  case  where  N  •+•  <*>  ,  that  is,  for  a  large  system  of 


equations.  The  eigenvalues  for  a  system  of  40  equations  are  verv 
close  to  the  limiting  vaLues.  The  eigenvalues  of  such  a  system 
are  bounded  above  bv  one.  The  minimum  eigenvalue  then  will 
control  the  st  :bi! Ity  of  the  system  and  is  given  when  n  =  N  by 


(xc) 

a 


l  +  lZpcc 

nfpTHj 


(13) 


Since  a  gives  the  "degree  of  implicitness"  of  equation  ill),  the 
entire  family  of  finite  element  expressions  can  be  investigated. 
Figure  1  shows  graphically  the  results  of  the  analysis  for  the 
pure- imp  L  ic  i  t ,  optimum-impLicit ,  Cranlc-Nicolson ,  and  explicit 
(Euler)  formulations.  It  is  immediately  obvious  that  the 
optimum-implicit  method  is  more  stable  than  the  Crank-Nicolson 
method.  This  is  not  the  case  with  finite  differences. 


RESULTS 

Extensive  calculations  were  made  using  equations  (11)  and 
(15).  Due  to  the  discontinuity  in  the  boundary  condition  at 
t  =  0  and  x  =  0,  problems  with  respect  to  convergence  of  the 
numerical  solution  were  encountered.  This  was  true  )lso  with  the 
Crandall  me:  h  d.  however.  To  eliminate  the  problems  owing  : 
the  discontinuity,  the  exact  analvtical  solution  to  equation  _ 
was  substituted  for  the  numerical  solution  at  the  end  of  the 
first  time  step.  The  convergence  of  the  numerical  solution  after 
a  change  in  the  discretization  scheme  from  10  to  20  elements  is 
shown  in  Figure  2.  The  discretization  error  ratio  (DER)  is 
defined  as  the  total  error  between  the  exact  and  numerical 
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solutions  for  E  elements  divided  by  the  error  at  the  same  point 
when  E  has  been  doubled.  Here  the  truncation  error  has  been 
summed  over  all  9  internal  nodes  for  the  10  element 
discretization  scheme  and  the  corresponding  9  internal  nodes  for 
the  20  element  scheme  and  is  termed  "generalized  mean  error."  In 
Figure  2,  the  DER  is  graphed  as  a  function  of  time  for  three 
values  of  the  parameter  a,  while  in  Figure  3,  it  is  graphed  as  a 
function  of  a  for  three  points  in  time.  The  elemental  truncation 
error  for  the  Crank-Nicolson  method  (a  =  .5)  and  the 
pure-implicit  method  (a  =  1.0)  is  approximately  1/4  if  Ax  is 
halved.  For  the  optimum  formulation  given  here,  the  elemental 
truncation  error  is  approximately  1/16  if  Ax  is  halved.  It  is 
clear  from  Figure  2  that  all  three  methods  approach  their 
theoretical  improvement  ratios.  It  is  also  interesting  to  note 
that  the  H°  error  norm,  defined  by 

II*  '  *EH0  =  /  4  ‘  *E)2  d*  (19) 

shows  a  peak  structure  in  Figure  4  similar  to  that  observed  for 
the  truncation  error  shown  in  Figure  3.  This  observation  leads 
to  speculation  that  there  may  be  a  predictable  value  of  a  which 
optimizes  the  finite  element  system  (11)  with  respect  to  the  H° 
error  norm. 

CONCLUSIONS 

1.  The  optimum  method  derived  here  using  Linear  finite  elements 
yields  the  same  system  of  equations  as  the  Crandall  method  for 
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pure  finite  differences.  The  finite  element  solution  is 
continuous,  however,  rather  than  discrete. 

2.  The  stability  for  the  optimum-implicit  scheme  is  greater  than 
for  the  Crank-Nicolson  method  since  the  oscillation  limits  are 
less  restrictive  than  in  the  latter  method. 

3.  The  theoretical  increase  in  accuracy  for  the  optimum  method 
over  the  Crank-Nicolson  method  is  clearly  verified  by  the 
numerical  results. 

4.  The  discontinuity  in  the  boundary  condition  at  t  =  0,  x  =  0 
for  the  Neumann  problem  causes  a  loss  in  the  order  of  convergence 
which  may  be  recoverable  by  using  a  discontinuous  elemental  trial 
function  at  the  point  of  the  discontinuity. 
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Figure  1.  Stability  Curves  for  Finite  Elements. 
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DISC.  ERROR  RATIO 

R).0  2.5  S.O  7.5  10.0  12.5 


®  Crank-Nicolson  Method,  a  =  0.5 
A  Optimum  Implicit  Method,  a  =  0.5417 

+  Pure  Implicit  Method,  a  =  I  .  0 

Fourier  Modulus:  p  =  2.0 

Discretization  Change:  Ax  =  0.1  to  ax  =  0.05 

Figure  2.  Convergence  as  a  Function  of  Time  for  the 
truncation  error. 
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Fourier  Modulus:  p  =  2.0 

Discretization  Change:  Ax  =  0.1  to  Ax  =  0.0S 
Figure  3.  Convergence  as  a  Function  of  a  in  the  L-,  Norm. 
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DISC.  ERROR  RflTIO 
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Fourier  Modulus:  p  =  2.0 

Discretization  Change:  =  0.1  to  Ax  =  O.OS 

Figure  4.  Convergence  as  a  function  of  a  for  the 
truncation  error. 
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